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1.  SUMMARY 


The  differences  between  earthquakes  and  explosions  are  largest  in  the  highest 
recordable  frequency  band.  In  this  band,  scattering  of  elastic  energy  by  small-scale 
heterogeneity  (less  than  a  wavelength)  can  equilibrate  energy  on  components  of  motion 
and  stabilize  the  behavior  of  the  Lg  wave  trapped  in  the  Earth's  crust.  Larger-scale 
detenninistic  structure  (greater  than  a  wavelength)  can  still  assume  major  control  over  the 
efficiency  or  blockage  of  the  Lg  and  other  regional/local  seismic  waves.  This  project 
models  the  combined  effects  of  the  large-scale  (detenninistic)  and  the  small  scale 
(statistical)  structure  to  invert  for  improved  structural  models  and  to  evaluate  the 
perfonnance  of  yield  estimators  and  discriminants  at  selected  International  Monitoring 
System  (IMS)  stations  in  Eurasia.  These  tasks  are  approached  by  synthesizing 
seismograms  using  a  radiative  transport  technique  to  predict  the  high  frequency  coda 
(>5  Hz)  of  regional  seismic  phases  at  stations  having  known  large-scale  three- 
dimensional  structure,  combined  with  experiments  to  estimate  the  effects  of  multiple¬ 
scattering  from  unknown  small-scale  structure.  We  describe  a  radiative  transport  code 
and  report  the  results  of  its  use  modeling  regional  wave  propagation  in  detenninistic  and 
statistical  3-D  structure,  including  the  computation  of  quantities  required  to  synthesize 
high  frequency  body  wave  coda  generated  by  small-scale,  statistically  described 
heterogeneity. 

2.  INTRODUCTION 

This  section  provides  some  background  and  introductory  material.  It  is  largely 
unchanged  from  our  earlier  report,  but  is  included  to  eliminate  the  need  for  readers  to 
refer  to  our  earlier  report.  In  Section  3  we  present  the  methods  and  procedures  used  to 
develop  the  code  with  updated  tests  and  evaluations.  Section  4  contains  test  our  recent 
results  and  discussion.  Concluding  remarks  are  presented  in  Section  5. 

2.1.  The  importance  of  very  high  frequency  coda  modeling 

Differences  in  source  properties  between  a  theoretically  ideal  earthquake  and 
explosion  are  well  understood:  distributed  slip  on  a  fault  plane  versus  an  approximate 
step-function  in  pressure  applied  to  a  spherical  cavity  deep  enough  for  explosion 
containment  (Brune,  1970;  Mueller  and  Murphy,  1971).  These  simple  theories  do  not 
incorporate  dynamically  triggered  fault  motion,  complex  surfaces  of  equivalent  explosion 
volumes,  and  explosive  tectonic  release.  Nevertheless  they  are  good  enough  to  accurately 
predict  that  the  scalar  moment  of  an  earthquake  will  have  much  higher  comer  frequency 
in  its  displacement  spectrum.  The  success  of  common  seismic  discriminants  based  on 
amplitude  ratios  in  different  group  velocity  windows  (e.g.,  P/Lg)  depend  as  much  on  this 
spectral  difference  as  they  do  on  differences  in  source  depth  and  P  and  S  radiation 
patterns  (Walter  et  al,  2009). 

High  frequency  Lg  coda  has  also  been  shown  to  be  useful  for  yield  estimation  (e.g., 
Mayeda,  1993).  As  frequency  increases,  however,  it  becomes  increasingly  difficult  to 
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obtain  perfect  knowledge  of  small-scale  structure  needed  to  predict  a  seismogram  of  an 
earthquake  or  an  explosion  wiggle  for  wiggle.  Yet  even  without  this  knowledge,  simple 
statistical  descriptions  of  medium  properties  can  be  used  to  predict  the  shape  of  high 
frequency  Lg  coda. 

2.2.  Limitations  of  numerical  modeling:  high  frequency  and  long  range 

Given  the  importance  of  high  frequencies  to  seismic  monitoring,  it  would  be  desirable 
to  model  high  frequency  regional  seismograms  with  enough  computational  efficiency  that 
many  structural  and  source  models  could  be  routinely  tested  and  compared  against 
observations  within  several  hours  or  less.  Unfortunately,  even  with  the  advance  of 
parallel  computation  and  increased  access  to  large  clusters  of  processors,  numerical 
modeling  in  three-dimensional  earth  models  rarely  exceeds  ranges  of  several  hundred 
wavelengths  for  computing  times  less  than  1  day.  At  5  Hz  this  range  is  typically  less  than 
100-200  km,  which  is  valuable  to  predictions  of  earthquake  strong  ground  motion,  but 
more  limited  in  value  to  nuclear  monitoring  where  IMS  stations  may  be  1000  km  or 
further  from  regions  of  interest  in  Eurasia. 

2.3.  Importance  of  large-scale  structure  and  limitations  of  a  calibration 
approach 

Given  the  apparent  success  of  discriminants  and  yield  estimators  from  high  frequency 
coda  with  minimal  knowledge  of  earth  structure  and  source  descriptions,  one  may  believe 
that  monitoring  is  purely  a  problem  of  calibrating  coda  measurements  and  amplitude 
ratios  along  paths  to  each  monitoring  station,  applying  spatial  interpolators  (e.g,  kriging) 
and  a  model  of  propagation  in  a  plane-layered  medium  (e.g.,  MDAC).  Many  examples 
have  been  collected  in  which  sharp  gradients  in  the  Moho  depth,  mountain  roots,  and 
deep  sedimentary  basins  have  strongly  affected  azimuthal  variations  in  the  character  and 
detection  of  regional  seismic  phases. 

The  effects  of  larger-scale  (greater  than  a  wavelength)  structure  will  limit  the 
portability  of  discriminants  and  yield  estimators  and  limit  the  effectiveness  of  their  path 
interpolations  and  calibrations.  In  addition  to  these  operational  concerns,  a  number  of 
fundamental  problems  involving  the  interaction  of  large  with  small-scale  structure  have 
not  been  well  investigated: 

•  Can  large-scale  structures  always  block  the  highest  detectable  frequency  of 
some  regional  phases? 

•  Is  there  a  frequency  band  on  Earth  in  which  the  effects  on  regional 
seismograms  of  large-scale  structure  disappear? 

•  How  important  is  small-scale  structure  in  the  vicinity  of  large-scale  structure? 

•  Is  anisotropy  of  the  scale-lengths  of  small-scale  structure  (horizontally  or 
vertically  stretched  heterogeneity)  important  to  regional  phase  propagation? 

By  combining  the  effects  of  large-scale  and  small-scale  structure  in  a  computationally 
efficient  3-D  method  of  seismogram  synthesis  based  on  the  radiative  transport  approach, 
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our  project  seeks  to  determine  small-and  large-scale  3-D  structure  in  the  vicinity  of 
selected  Eurasian  IMS  stations  by  matching  the  coda  of  very  high  frequency  (>  5  Hz) 
regional  phases.  The  results  of  this  modeling  will  be  used  to  optimize  coda-based  yield 
estimators  and  discriminants  as  a  function  of  azimuth  to  these  stations. 

3.  TECHNICAL  APPROACH 

3.1.  Radiative  Transport 

Radiative  transport  seeks  to  model  the  diffusive  transport  of  energy  through  processes 
of  multiple  scattering  in  a  medium  having  random  fluctuations  in  physical  properties.  Its 
earliest  and  most  extensive  applications  are  to  the  propagation  of  electromagnetic  waves, 
especially  to  the  description  of  stellar  images  and  scintillation  (Chandrasekhar,  1960). 

Wu  (1985)  introduced  radiative  transport  to  seismology,  modeling  the  coda,  or  envelope 
of  energy,  starting  and  following  elastic  P  and  S  waves.  By  including  multiple-scattering, 
radiative  transport  improved  upon  the  single  scattering,  Rayleigh-Bom,  model  of  seismic 
coda  introduced  by  Aki  and  Chouet  (1975).  An  early  significant  result  of  Wu’s  work  with 
radiative  transport  is  that  multiple  scattering  often  becomes  the  dominant  mechanism 
controlling  coda  shapes  at  frequencies  above  5  Hz  on  Earth,  explaining  a  spindle-shaped 
(lunar-like)  coda  observed  in  this  frequency  band.  Reviews  of  seismic  radiative  transport 
modeling  include  the  text  by  Sato  and  Fehler  (1998)  and  a  monograph  by  Margerin 
(2004).  Przybilla,  et  al.  (2009)  has  compared  the  method  against  finite-difference 
synthetics  in  2-D,  validating  body  wave  coda  predictions  of  the  radiative  transport 
method. 

3.2.  Construction  of  models:  statistically  described  structure 

Frankel  and  Clayton  (1986)  demonstrated  how  earth  models  can  be  constructed  that 
have  a  specified  heterogeneity  spectrum.  The  procedure  consists  in  having  perturbations 
at  knot  points  driven  by  a  random  number  generator,  Fourier  transfonning  the  spatial 
model  into  the  wavenumber  domain,  filtering  by  a  desired  wavenumber  spectral  shape, 
and  inverse  Fourier  transfonning  back  to  space. 

Since  Frankel  and  Clayton’s  original  work,  advances  have  occurred  in  statistical 
characterization  and  understanding  of  small-scale  geologic  heterogeneity.  Goff  et  al. 
(1994)  described  a  procedure  by  which  models  having  the  spatial  statistics  of 
polycrystalline  or  multi-modal  assemblages  of  rocks  can  be  generated.  Work  beginning 
with  Levander  et  al.  (1994)  and  Pullammanappallil  et  al.  (1997)  makes  it  possible  to 
fonnulate  statistical  models  for  common  sedimentary  and  metamorphic  formations. 
Model  statistics  can  be  modified  as  needed  to  reproduce  the  three-component  behavior  of 
the  coda  of  regional  seismic  phases.  Strong  emphasis  is  placed  on  correctly 
characterizing  the  statistics  of  small-scale  heterogeneities  in  the  upper,  highly 
heterogeneous,  5  km  of  the  earth.  Small-scale  statistics  affect  the  partitioning  of  P  and  S 
energy  on  the  three-components  of  motion  by  the  effects  of  scattering  near  the  source  and 
receiver  and  hence  the  perfonnance  of  discriminants  is  tuned  to  differences  in  sources 
occurring  at  the  depths  of  contained  nuclear  tests. 
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3.3.  Construction  of  models:  deterministic  structure 

Lateral  variations  in  crustal  thickness,  basin  depths,  mountain  roots,  and  lateral 
tectonic  transitions  significantly  affect  the  phases  used  for  discrimination  and  detection. 
One  example  is  the  study  by  Pedersen  et  al.  (1998),  who  explain  anomalous  Rayleigh  to 
Love  mode  conversion  from  Lop  Nor  explosions  by  changes  in  crustal  thickness  at  the 
boundary  of  the  Tarim  Basin  and  Tian  Shan  mountain  belt.  Moho  topography  and  basin 
thickness  can  also  strongly  affect  the  propagation  of  Lg  (Connier  and  Anderson,  2004). 
The  scale  of  these  types  of  lateral  structural  variations  is  often  large  enough  to  be 
resolvable  by  local  and  regional  reflection  and  refraction  experiments,  gravity  and 
magnetic  data,  regionalization  by  surficial  geology,  and  global  surface  wave  inversions. 
Hence,  we  refer  to  these  types  of  structures  as  detenninistic.  The  types  of  data  used  to 
infer  detenninistic  structure  are  collected  at  widely  different  spatial  scales,  presenting  a 
challenge  to  the  parameterization  of  a  three-dimensional  model  appropriate  for  a  region 
surrounding  a  particular  seismic  station.  The  parameterization  should  be  flexible  enough 
to  be  specified  at  high  resolution  where  data  is  available  and  at  lower  resolution  where  it 
is  not.  Resolution  should  be  high  to  describe  features  important  to  regional  wave 
propagation,  such  as  Moho  and  basin  topography,  but  can  be  lower  near  interfaces  having 
smaller  velocity  contrasts  and  lower  with  increasing  depth  in  the  mantle,  where 
heterogeneity  power  decreases. 


3.4.  Synthetic  Algorithm 

We  are  developing  a  3-D  radiative  transport  code  for  synthetic  seismogram 
generation  called  Radiative3D.  Radiative3D  is  a  new  code  that  builds  upon  and 
augments  two  existing  codes,  namely  Raytrace3d  (Menke,  2002)  and  PSPhonon  (Shearer 
and  Earle,  2004;  2008).  Raytrace3D  solves  ray  propagation  in  a  3-D  velocity  model,  and 
PSPhonon  performs  radiative  transport,  including  scattering,  in  a  1-D  layered-earth 
model.  Our  code  performs  both  tasks  in  a  full  3-D  model. 

Radiative  transport  theory  treats  seismic  energy  as  composed  of  discrete  energy 
packets  leaving  elements  of  solid  angle  surrounding  a  source  region,  and  uses  ray  theory 
to  trace  the  paths  of  these  energy  packets  (phonons)  through  a  detenninistic  velocity 
model.  Travel  time,  path  integrated  intrinsic  attenuation,  and  polarization  of  particle 
motion  are  tracked  as  the  phonons  propagate.  Seismic  energy  is  finally  collected  and 
tabulated  into  temporal  and  spatial  bins  when  they  reach  the  surface  of  the  model  or  some 
other  point-of-interest  in  the  model.  This  allows  the  generation  of  synthetic  seismograms 
at  specified  model  locations. 

The  radiative  transport  approach  has  several  advantages  and  some  short-comings 
which  differentiate  it  from  other  approaches.  In  particular,  radiative  transport: 

•  Requires  “spraying”  a  sufficient  quantity  of  phonons  from  the  source  event, 
covering  a  range  of  take-off  angles  designed  to  mimic  either  earthquake  or 
explosion  energy  release. 
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•  Is  generally  faster  than  solving  the  full  wave  equation. 

•  Uses  ray  theory  to  model  energy  transport  through  large-scale  “deterministic” 
structure  (valid  when  wavelengths  are  sufficiently  small  compared  to  length 
scales  of  the  deterministic  structure  being  modeled). 

•  Effects  of  small-scale  structure  are  modeled  using  Monte  Carlo  methods  by 
implementing  a  statistical  scattering  model. 

•  The  Lg  wave  and  its  coda  are  modeled  by  S  body  waves,  multiply  reflected  at 
critical  incidence  from  the  Moho,  multiply  scattered  and  trapped  in  the  crust. 

•  The  Pn  wave  and  its  coda  are  modeled  by  P  waves  interacting  with  and  turning 
beneath  the  Moho  and  scattered  by  statistical  structure  surrounding  the  Moho. 

•  The  Rg  wave  cannot  be  properly  represented  as  a  fundamental  model  mode 
Rayeligh  wave  having  frequency-dependent,  exponentially  decaying,  energy 
beneath  the  free  surface.  This  limitation,  however,  is  mitigated  by  observations 
that  very  high  frequency  (>2  Hz)  Rg  is  normally  absent  in  local/regional 
seismograms  beyond  200  km  range  due  to  strong  scattering  by  surface  topography 
along  its  path. 

The  development,  testing,  and  validation  of  Radiative3D  includes  a  number  of  steps, 
which  are  discussed  below  with  their  related  results.  Steps  discussed  in  last  year’s  report 
include: 

•  The  ability  to  spray  P-  and  S-polarized  phonons  from  a  defined  source  according 
to  a  pattern  specified  by  a  six-element  moment  tensor,  which  can  describe  both 
explosion  and  shear-dislocation  point-source  events; 

•  The  ability  to  propagate  phonons  through  a  whole-space  elastic  medium  of 
constant  seismic  velocity; 

•  The  ability  to  simulate  a  scattering  event  after  propagating  a  randomly  chosen 
distance  selected  from  a  mean-free-path  distribution  detennined  by  the  scattering 
parameters  of  the  medium; 

•  The  ability  to  modify  the  trajectory  and  polarization  of  the  phonon  at  a  scattering 
event  according  to  the  scattering  functions  of  Sato  and  Fehler  (1998); 

•  The  ability  to  collect  and  record  the  seismic  energy  of  the  phonons  at  plane 
interfaces  located  a  specified  distance  from  the  source  and  resolve  it  into  three 
components  of  motion,  allowing  the  generation  of  synthetic  seismic  envelopes. 

3.4.1.  Visualizing  Output 

Radiative3D  produces  output  suitable  for  visualization  with  external  tools.  The  output 
capabilities  fall  into  two  major  categories:  (1)  event  reporting,  and  (2)  seismic  energy 
binning.  Event  reporting  means  reporting  the  progress  of  individual  phonons  in  a  play- 
by-play  manner.  Examples  of  “events”  in  this  context  include:  generation,  scattering, 
crossing  a  model  boundary,  free-surface  or  discontinuity  interface  reflection  and 
transmission,  etc.  When  these  events  are  analyzed  in  post-processing,  detailed  pictures  of 
energy  propagation  can  be  constructed.  One  of  the  first  visualization  tools  developed  was 
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a  GNU  Octave  script  to  produce  videos  illustrating  P  and  S  energy  propagation  in  the 
Earth  model. 


Figure  1  shows  energy  propagation  visualized  as  a  time  series  plot  of  phonon  locations  in 
a  model  cross  section.  Red  dots  represent  P  phonons,  blue  dots  represent  S  phonons. 
Reflection  and  refraction  occurs  at  interfaces  between  model  layers  at  which  velocities 
are  discontinuous.  Conversions  between  P  and  S  polarization  are  driven  by  both 
scattering  and  reflection/transmission  events. 


Earthquake  Time-Series: 


Explosion  Time-Series: 


Al  T=7-° 


Figure  1.  Energy  propagation  visualized  as  a  time  series  plot  of  phonon  locations  in  a 
model  cross  section.  Note  the  intial  absence  of  the  S  wavefron  ts  (blue)  in  the 
explosion  time-series  and  its  progressive  development  with  time  due  to  multiple 
scattering  and  interaction  with  the  free  surface. 


3.4.2.  Seismometer  Output 

The  other  fonn  of  output  which  can  be  analyzed  for  quantification  or  visualization  is 
binned  seismic  energy.  Radiative3D  allows  for  the  placement  of  virtual  seismometers 
along  any  surface  designated  as  a  collection  surface.  Whenever  a  phonon  interacts  with  a 
collection  surface  within  a  specified  gather  radius  of  a  seismometer,  the  energy  of  that 
phonon  is  collected  by  the  seismometer,  decomposed  into  three  cartesian  energy  axes, 
and  binned  into  time  windows.  At  the  end  of  simulation,  these  energy  bins  are  output  and 
can  be  interpreted  as  seismic  energy  traces,  suitable  for  making  envelope  plots.  Traces 
from  multiple  seismometers  arranged  in  a  linear  array  can  be  combined  to  make  travel¬ 
time  curves. 
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Figure  2  shows  two  seismic  envelope  plots  produced  by  the  Radiative3D  code.  The 
traces  represent  4.0  Hz  phonon  energy  arriving  at  a  virtual  seismometer  approximately 
900  km  from  a  hypothesized  earthquake  source  (left)  and  explosion  source  (right) 
synthesized  in  a  3-D  model  of  the  Lop  Nor  region  having  a  simple  deterministic  and  trial 
statistical  structure. 


Location:  887.61  km  at  276.24  deg  az  from  source 


300 

Time  (s) 


Figure  2.  Seismic  envelope  plots  produced  by  the  Radiative3D  code.  Traces  are  shown 
for  components  of  motion  NS  (Y),  EW  (X),  and  vertical  (Z)  and  wave  type  (P  or  S) 
arriving  at  the  receiver. 


3.4.3.  New  Features  Complete: 

In  addition  to  features  reported  in  last  year’s  report,  the  following  new  features  are  fully 
or  substantially  complete  as  of  this  writing: 

Layered  multi-cell  Earth  models: 

While  full-3D  Earth  models  based  on  tetrahedral  model  cells  are  not  yet  supported,  an 
intermediate  cell  type  based  on  stacked  cylinders  is  supported  and  functional  for 
simulating  layered-Earth  models.  Partial-3D  functionality  is  available  because  the 
interfaces  between  adjacent  cylinders  can  be  given  an  arbitrary  inclination.  This  means 
that  first-order  lateral  variations,  such  as  a  sloped  Moho  can  be  modeled. 

Reflection/Transmission  scattering  at  discontinuity  interfaces: 

Free-surface  reflections,  as  well  as  sharp  velocity  discontinuities  (such  as  at  the  Moho 
layer),  require  special  treatment  of  ray  propagation.  Much  like  the  radiation  patterns  of 
scattering  from  small-scale  volumetric  inhomogeneities,  a  ray’s  interaction  with  an 
interface  across  which  seismic  velocities  and/or  material  densities  abruptly  change  can  be 
thought  of  as  a  form  of  scattering  radiation  pattern.  For  a  given  incident  ray,  there  are  up 
to  four  different  outgoing  ray  directions,  and  two  different  ray  types  (P  and  S),  with 
special  care  being  necessary  to  correctly  transform  the  S-polarization  angle.  Aki  and 
Richards  (1980)  develop  a  set  of  20  scattering  coefficients  that  cover  all  possibilities 
between  incoming  P  and  S  waves  and  their  associated  transmitted  or  reflected  P  or  S 
outgoing  waves.  These  coefficients  quantify  the  amplitude  and  phase  of  the  outgoing 
waves,  which  is  used  by  Radiative3D  to  compute  probabilities  of  an  incoming  ray 
transforming  into  the  corresponding  outgoing  ray  type.  This  allows  for  the  modeling  of 
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reflection  (at  a  free  surface  or  interface)  and  transmission  (interface)  of  rays  with 
corresponding  possibility  of  transfonnation  between  P  and  S  polarization  types.  This 
feature  has  been  recently  completed  in  the  Radiative3D  codebase. 

3.4.4.  Some  Notable  Features  Not-Yet  Complete: 

Software  development  is  ongoing,  and  while  we  have  reached  a  state  where  certain 
studies  of  real-world  significance  can  be  undertaken  with  the  code  as-is,  there  still  remain 
some  features  to  be  developed  before  we  reach  full  target  functionality.  The  most 
significant  of  these  are: 

Tetrahedral  Model  Cells: 

Three-dimensional  Earth  models,  having  fully  general  spatial  gradients  of  velocity  and 
density  within  layers  need  to  be  discretized  in  some  fashion.  The  method  we  have  chosen 
is  to  divide  the  model  into  tetrahedral  cells  inside  of  which  the  Earth  parameters  are 
analytically  specified.  A  tetrahedron  is  the  polyhedron  with  the  minimum  number  of 
vertices  needed  to  specifiy  a  three-dimensional  volume,  which  allows  for  straighforward 
approaches  to  dividing  space  into  cellular  components.  Another  advantage  of  tetrahedra, 
(over,  e.g.,  rectangular  prisms),  is  that  Earth  parameters  specified  on  the  four  vertices 
allows  unique  and  complete  specification  of  the  properties  in  terms  of  linear  gradients 
throughout  the  volume  of  the  cell. 

Much  of  the  infrastructure  for  cellular  model  discretization  is  already  complete  in  the 
code,  and  is  utilized  in  the  layer-celled  structure  already  supported  in  the  code. 
Implementing  cells  with  tetrahedral  geometry  remains  primarily  an  implementation  detail 
at  this  point. 

Velocity  Gradients: 

In  a  unifonn-velocity  medium,  ray  paths  are  straight  lines.  In  the  real  Earth,  rays  can 
follow  curvilinear  paths.  The  plan  for  Radiative3D  is  for  velocities  to  be  specified  as 
first-order  gradients,  resulting  in  paths  which  are  circular  arcs  withing  a  given  model  cell. 
This  allows  for  a  much  more  natural  representation  of  real  Earth  structure,  and  lessens  the 
detail  needed  to  represent  that  structure.  At  present  Radiative3D  supports  only  unifonn 
velocities  (gradients  on  the  macro  level  can  be  represented  as  a  multicellular  stair-step 
model).  The  next  development  step  of  Radiative3D  will  be  to  implement  support  for 
velocity  gradients.  This  step  will  naturally  follow  the  implementation  of  tetrahedral 
model  cells,  since  this  geometry  allows  the  unique  specification  of  the  gradients  by 
specifying  known  velocities  on  the  cell  vertices. 

3.4.5.  Features  Not  Currently  Targeted: 

In  addition  to  the  code  features  detailed  above,  another  feature,  not  currently  targetted  for 
development,  remains  within  feasible  reach  of  the  codebase  as  it  now  stands,  and  may 
provide  a  starting  point  for  future  work.  In  particular,  while  Radiative3D  currently 
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produces  seismic  envelope  plots,  it  is  technically  possible  to  produce  full-waveform 
seismograms  by  tracking  not  only  the  amplitude,  but  also  the  phase,  of  each  phonon  as  it 
propagates  through  the  model.  By  binning  the  seismic  energy  as  amplitude  phasors 
rather  than  energy  scalars,  full  seismograms,  rather  than  envelopes,  can  be  produced. 
Phase  tracking  through  direct  propagation  is  trivial,  as  it  is  proportional  to  the  product  of 
frequency  and  travel  time.  Slightly  more  intricate  are  the  phase  shifts  that  result  from 
scattering  events,  and  from  reflection/transmission  interactions  with  discontinuity  or  free- 
surface  interfaces.  However,  these  phase  factors  are  actually  already  calculated  by  the 
codebase  when  it  calculates  scattering  and  reflection/transmission  probabilities.  Adding 
the  accounting  necessary  to  track  cumulative  phase  through  a  phonon  lifetime  would  be  a 
very  tractable  problem. 

4.  RESULTS  AND  DISCUSSION 
4.1.  Earth  Model  for  Region  of  Interest 

We  originally  planned  to  focus  our  study  on  the  area  in  the  vicinity  of  the  NIL 
seismic  station.  We  discovered  that  while  there  were  many  earthquakes  with  regional 
seismograms  in  that  area,  regional  explosion  data  were  severely  lacking.  We  next 
planned  to  use  regional  explosion  data  from  the  Semipalatinsk  Test  site  (STS)  in 
Kazakhstan  as  received  at  the  Borovoye  (BRV)  Geophysical  Observatory  in  Kazakhstan. 
Currently  there  are  few  earthquake  waveforms  in  that  area  to  compare  with  explosion 
waveforms.  (This  may  be  alleviated  in  the  near  future  by  the  LDEO  digitization  project 
of  analog  seismograms  recorded  by  regional  and  PNE  associated  networks  in  the  former 
USSR).  We  decided  to  switch  to  the  Lop  Nor  explosion  area  since  there  were  both 
explosion  and  earthquake  waveforms  in  the  same  area,  in  addition  to  the  availability  of 
detailed  maps  of  regional  wave  propagation  efficiencies  by  LANL.  Lop  Nor  is  located 
near  the  southeastern  side  of  the  Tien  Shan,  a  region  of  moderate  earthquake  activity  and 
contemporary  horizontal  compressive  stress  in  the  earth’s  crust,  which  proves  to  be  of 
significance  in  the  discrimination  of  earthquakes  from  nuclear  explosions  in  that  region. 
Sykes  and  Nettles  (2009)  found  that  more  than  half  of  the  total  numbers  of  earthquakes  in 
the  Reviewed  Event  Bulletin  (REB)  of  the  IMS  that  occurred  within  100  km  (62  mi)  of 
six  test  sites  from  2000  through  2008  occurred  near  Lop  Nor. 

Figure  3  shows  the  area  we’re  currently  investigating.  The  Chinese  station  Urumqi 
(assigned  the  station  code  WMQ  by  seismologists)  is  about  2.15°  (250  km  or  156  mi) 
from  Lop  Nor.  MAK  and  WUS  are  ~6.85°  from  Lop  Nor.  For  the  earthquake  and 
explosion  data  we  use  seismograms  from  events  found  using  the  Incorporated  Research 
Institutions  for  Seismology  (IRIS)  web  site. 
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Figure  3.  Lop  Nor  area  and  some  regional  seismic  stations. 


Although  we  have  not  constructed  a  detailed  3-D  earth  model  for  this  region,  we 
constructed  a  simple  3-D  model  to  begin  synthetic  tests.  The  model  was  constructed  of 
layers  of  uniform  velocity,  with  dipping  planes  separating  the  layers  (Figure  4).  Current 
functionality  in  Radiative3D  allows  these  interface  planes  to  take  on  arbitrary  3-D 
orientation.  Depth  profiles  from  CRUST2.0  topographic  elevations,  and  Moho  depths  at 
three  locations  (Lop  Nor,  MAK,  and  WUS)  were  used  to  locate  and  orient  the  planes. 


LOP 


MAK 


LOP 

tSits 


WUS 


IF 


Upper  Crust 

Scattering  Scale  Length: 
a  =  0.50  km 

Middle  Crust 

Scattering  Scale  Length: 
a  =  0.50  km 

Lower  Crust 

Scattering  Scale  Length: 
a  =  0.50  km 

Mantle 

Scattering  Scale  Length: 
a  =  1.00  km 

Figure  4.  Simplified  Lop  Nor  earth  model.  Note  a  thin  sediment,  low  P  and  S  velocity, 
layer  is  included. 
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4.2.  Regional  Data 

Figure  5  shows  some  waveform  examples  from  the  Lop  Nor  region.  The  data  come  from 
an  explosion  at  the  Lop  Nor  site  as  recorded  at  station  MAK  at  great  circle  distance  6.85° 
and  an  earthquake  in  the  same  area  recorded  at  both  MAK  and  WUS.  WUS  is 
approximately  at  the  same  distance  from  Lop  Nor  as  MAK  but  in  at  different  azimuth. 
The  left  column  in  Figure  5  shows  the  data  bandpass  filtered  1-2  Hz  and  the  right  column 
shows  filtering  6-8  Hz.  Not  only  are  there  differences  between  the  earthquake  and 
explosion  data  at  MAK,  there  are  also  differences  between  the  earthquake  traces  at  MAK 
and  WUS.  We  plan  to  use  our  synthetics  to  explore  and  understand  these  differences 
from  effects  of  both  large-  and  small-scale  3-D  structure. 


1-2  Hz.  Explosion.  Lop  Nor-MAK 


SoUtwrUlryan 
MAK  BHZ 
OCT  Of  <?»).  It* 
032(42017 
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6-8  Hz.  Explosion.  Lop  Nor-MAK 


1-2  Hz.  Earthquake.  Lop  Nor-MAK 


6-8  Hz,  Earthquake,  Lop  Nor-MAK 


1  -2  Hz,  Earthquake,  Lop  Nor-WUS 


6-8  Hz,  Earthquake,  Lop  Nor-WUS 


Figure  5.  Lop  Nor  regional  vertical  component  bandpassed  waveforms. 


4.3  Synthetic  Results 

We  have  produced  envelope  and  travel-time  synthetics  in  Radiative3D  using  a 
simplified  Earth  model  serving  as  a  stand-in  for  more  complex  models  to  be  supported  in 
future  revisions  as  work  on  Radiative3D  comes  to  completion.  The  Earth  model  used  is  a 
layered  model  consisting  of  layers  of  spatially  unifonn  parameters  (velocity,  density, 
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scattering  parameters).  The  planar  interfaces  separating  the  layers  need  not  be  parallel, 
and  in  fact  can  take  on  arbitrary  orientations.  We  have  taken  advantage  of  this  freedom  to 
construct  a  first-order  approximation  of  the  topography  (both  surface  and  subsurface, 
e.g.,  Moho)  in  the  region  of  the  Lop  Nor  nuclear  test  site.  Using  known  elevations  of 
Lop  Nor  and  seismic  stations  MAK  and  WUS,  along  with  Moho  depths  from  the  Cornell 
Moho  model  and  layer  profiles  from  CRUST2.0  at  those  same  locations,  we  located  and 
oriented  a  set  of  four  crustal  layers  (sediments,  upper,  middle,  and  lower  crust,  and  the 
top  layer  of  the  mantle).  We  also  defined  an  additional  16  mantle  layers  from  AK-135,  to 
a  depth  of  859  km.  These  layers  served  as  a  model  in  which  to  run  early  test  synthetics. 

The  model  assumed  a  minimum  characteristic  scale  length  in  an  exponential 
autocorrelation  of  heterogeneities  of  0.50  km  in  the  crust,  1.00  km  in  the  mantle,  and 
0.25  km  in  the  thin  sediments  layer.  (An  exponential  autocorrelation  provides  a  white 
heterogeneity  spectrum  up  to  the  minimum  scale  length,  at  which  the  heterogeneity 
power  decays  as  k'2  with  increasing  wave  number  k).  Scattering  strengths,  defined  as  the 
fractional  velocity  perturbations  that  characterize  the  heterogeneity,  were  0.01  in  the 
crust,  0.008  in  the  mantle,  and  0.012  in  the  sediments  layer.  Figure  4  shows  cross- 
sections  illustrating  the  crustal  layers  of  the  test  model  along  the  two  paths. 

Situated  within  the  model  we  place  two  arrays  of  160  virtual  seismometers  each  along 
the  paths  of  Lop  Nor  to  WUS  and  Lop  Nor  to  MAK.  These  arrays  collect  and  bin  the 
phonon  counts  to  produce  synthetic  envelopes  and  travel  time  curves. 

Synthetic  envelopes  are  for  a  frequency  of  4.0  Hz  and  an  assumption  that  all  phonons 
originate  at  time  t  =  0.  (Temporally  complex  events  can  be  simulated  in  post-processing 
by  convolution  with  a  source  time  function.)  We  modeled  two  distinct  source 
mechanisms  parameterized  by  moment  tensor  elements.  One  mechanism  was  a  pure 
isotropic  event  meant  to  simulate  an  idealized  explosion,  and  the  other  was  a  pure  shear 
event  to  simulate  an  earthquake.  We  set  the  depth  of  the  explosion  event  to  2.0  km  below 
the  surface,  and  the  depth  of  the  earthquake  event  to  42.0  km  below  the  surface, 
comparable  to  some  of  the  larger  magnitude  earthquakes  to  recorded  in  the  region.  The 
simulation  consisted  of  the  spraying  of  200  million  phonons  from  the  event  source  at 
probabilistically  chosen  take-off  angles  selected  from  a  discrete  set  of  5.2  x  106  possible 
take-off  angles  roughly  evenly  spaced  throughout  a  unit  sphere  (no  preferred  direction). 
The  source  radiation  pattern  is  determined  via  relative  probability  amplitudes  for  each 
take-off  angle  computed  from  moment  tensor  elements  at  program  initialization. 

Figure  6  shows  synthetics  from  the  earthquake  along  the  two  paths.  The  travel-time 
plots  show  two  major  divisions  of  energy  arrival,  with  the  lowest  slope  (time/range) 
arrival  representing  the  Pn-wave  arrival  and  the  higher  slope  arrival  representing  the  Sn- 
wave  and  succeeding  Lg  coda.  We  believe  the  banding  in  each  arrival  is  the  result  of 
surface  to  Moho  multiples  comprising  Lg,  but  this  is  unconfirmed  as  yet.  There  is  a  stark 
difference  in  the  relative  amplitude  of  the  P  and  S  arrivals  between  the  two  paths.  We 
believe  the  majority  of  this  difference  is  explainable  via  the  source  radiation  pattern  and 
that  only  minor  visible  differences  are  explained  by  topographical  differences  between 
the  Lop  Nor  to  WUS  and  Lop  Nor  to  MAK  paths. 
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Location:  887.61  km  at  276.24  deg  az  from  source 


Travel-time  curve  at  27 6. 24 -degrees  azimuth 


0  100  200  300  400  500  600  700  800 

Range  (km) 


0.2 


■  0.1 
I  0.05 


Location:  914.24  km  at  319.43  deg  az  from  source 


10C  200  300  400  500  6C0  700  800  900 
Range  <«*) 


0.15 


.1 


Figure  6.  Synthetics  from  an  earthquake  event  to  MAK  and  WUS. 


Figure  7  shows  synthetics  from  the  explosion  event  along  the  same  two  paths.  The 
most  noticeable  difference  between  these  and  the  earthquake  synthetics  is  the  absence  of 
a  coherent  S-wave  front.  Also  absent  is  any  obvious  banding  in  the  P-wave  arrival, 
which  is  curious,  but  may  be  a  result  of  the  much  shallower  source  depth,  resulting  in  a 
stronger  direct  arrival  and  suppression  of  surface  multiples  due  to  oblique  incidence  of 
the  surface  reflections.  Differences  between  the  two  paths  are  harder  to  see  in  the  travel¬ 
time  plots,  but  can  be  seen  more  readily  in  the  envelope  plots.  Since  the  event  source  is 
isotropic,  any  differences  in  energy  arrival  between  the  paths  must  result  from  path 
differences  in  topography.  The  path  from  Lop  Nor  to  MAK  is  more  sloping  and  the  path 
from  Lop  Nor  to  WUS  is  more  flat. 
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Location:  887.61  km  at  276.24  deg  az  from  source 


Travel-time  curve  at  27 6. 24 -degrees  azimuth 


Location:  914.24  km  at  319.43  deg  az  from  source 


Travel-time  curve  at  319.432-degrees  azimuth 
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Figure  7.  Synthetic  seismograms  from  an  explosion  recorded  at  MAK  and  WUS. 


5.  CONCLUSIONS 

We  have  continued  development  of  a  code  to  shoot  body  wave  rays  through  general 
deterministic  3-D  structure,  including  the  computation  of  quantities  required  to 
synthesize  high  frequency  body  wave  coda  generated  by  small-scale,  statistically 
described  heterogeneity.  Modifications  in  this  period  include  layer  boundaries  having 
generalized  dips  and  the  treatment  of  reflection-transmission  at  3-D  boundaries.  Features 
to  be  added  in  the  next  reporting  period  will  include  generalized  3-D  spatial  gradients 
within  layers,  and  the  effect  of  path  integrated  intrinsic  attenuation.  Spatial  velocity 
gradients  beneath  the  Moho  will  be  particularly  important  for  predicting  the  efficiency  of 
Pn  and  Sn  propagation. 

Coda  envelopes  can  be  plotted  in  several  styles  with  different  emphases:  as  a  function 
of  time  and  component  of  motion  at  a  single  station,  as  a  function  of  distance  and  time  to 
assist  in  travel-time  picking  and  travel-time  uncertainty  estimates,  or  as  snapshots  in  time 
as  a  function  of  depth  and  range  to  understand  the  evolution  of  P  and  S  energy  and  the 
homogenization  of  the  radiation  pattern. 

The  availability  of  both  explosion  and  earthquake  waveforms,  detailed  maps  of 
regional  phase  efficiency,  and  the  existence  of  many  well-constrained  3-D  deterministic 
structural  models  derived  from  transportable  array  studies  makes  the  Lop  Nor  region  an 
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ideal  region  to  implement  our  modeling  technique.  Our  goals  will  include  the  modeling 
of  effects  of  the  unknown  small-scale  structure  and  its  associations  with  regional  geology 
and  tectonics,  and  an  effort  to  explain  the  LANL  mapped  patterns  of  regional  phase 
propagation  efficiency. 
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